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Abstract 


A Legendre pseudo spectral philosophy based multi-phase constrained fuel-optimal trajectory design approach is presented in this 
paper. The objective here is to find an optimal approach to successfully guide a lunar lander from perilune (18 km altitude) of a transfer 
orbit to a height of 100 m over a specific landing site. After attaining 100 m altitude, there is a mission critical re-targeting phase, which 
has very different objective (but is not critical for fuel optimization) and hence is not considered in this paper. The proposed approach 
takes into account various mission constraints in different phases from perilune to the landing site. These constraints include phase-1 
(‘braking with rough navigation’) from 18 km altitude to 7 km altitude where navigation accuracy is poor, phase-2 (‘attitude hold’) 
to hold the lander attitude for 35 sec for vision camera processing for obtaining navigation error, and phase-3 (‘braking with precise 
navigation’) from end of phase-2 to 100 m altitude over the landing site, where navigation accuracy is good (due to vision camera nav- 
igation inputs). At the end of phase-1, there are constraints on position and attitude. In Phase-2, the attitude must be held throughout. At 
the end of phase-3, the constraints include accuracy in position, velocity as well as attitude orientation. The proposed optimal trajectory 
technique satisfies the mission constraints in each phase and provides an overall fuel-minimizing guidance command history. 
© 2017 COSPAR. Published by Elsevier Ltd. All rights reserved. 
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1. Introduction can serve as a base for solar power harvesting (Hickman 


et al., 1990). Further, it can also act as a base for launching 


Because of its proximity to earth, moon is often favored 
as a base for conducting new demonstrations in space tech- 
nology. Moreover, exploration of moon has attracted the 
global attention since it is found to be mineral rich and also 
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ambitious missions to reach halo orbits around a Lagran- 
gian point and to extending human explorations to inter- 
planetary destinations (Dunham et al., 2013). In addition, 
proof of existence of water on moon confirmed by the 
impact probe of Chandrayaan-1 (Pieters et al., 2009) has 
given enormous hope for assistance when an envisaged 
lunar research base is constructed on the surface of the 
moon. Hence, there is a renewed interest across the globe 
for thorough exploration of moon. 

Even though human missions have been carried out in 
the past for technology demonstrations, a more meaningful 
prolonged and cost-effective exploration of moon can hap- 
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pen only when the mission is carried out in autonomous 
mode with a lander-rover combination. A critical difficulty 
of such an autonomous mission, however, is the fact that 
the lander carrying the rover must soft land on the moon 
surface with the help of on-board sensors, processors and 
actuators without any human intervention. In order to 
meet the mission goals and precise landing sequence, the 
lunar landing trajectory has built-in terminal constraints 
at various points which must be satisfied. Landing missions 
on celestial bodies without atmosphere (such as moon) also 
demand careful planning and execution of engine burns, 
which minimizes the fuel consumption. If not optimized, 
this extra fuel must be carried in the spacecraft leading to 
lesser useful payload mass, and therefore, it is imperative 
that such a mission must be executed with minimum fuel 
requirement. This planning makes the overall mission quite 
cost effective. 

Many researchers in the past have attempted to land a 
probe on the moon and their results have become basis 
for some of the decisions in descending trajectory design. 
It is opt to provide their research contributions: An implicit 
guidance logic for soft lunar landing has been attempted 
(McInnes and Radice, 1996) in which the lander tries to 
reach the designated landing point by riding on the line 
of sight and gradually reducing the approach velocity. 
However, such an intuitive and simplistic strategy normally 
does not lead to a proper orientation of the vehicle at the 
landing site. More importantly, it does not lead to a fuel 
optimal trajectory. Based on optimal control theory, a gen- 
eric explicit guidance law has been proposed using lin- 
earized dynamics that minimizes the final time as well as 
the total acceleration in all three directions, thereby leading 
to a fuel optimal solution (DSouza, 1997). Moreover, the 
terminal boundary conditions on position and velocity 
imposed as hard constraints ensure that these constraints 
are met in an efficient way. Additionally, it leads to a closed 
form analytic formula for the necessary acceleration in 
both lateral and axial directions, which leads to ease of 
onboard implementation. Because of these characteristics, 
the method has attracted many researchers across the globe 
to use this approach to address vertical orientation of the 
lander while touchdown. 

The desired lander orientation demands terminal accel- 
eration which has to be enforced as a hard constraint in 
all three directions. The idea of including the acceleration 
as one of the state for addressing this constraint has been 
analyzed by (Uchiyama et al., 2005; Mathavaraj and 
Padhi, 2017). Even though, the proposed analytical guid- 
ance logic ensured vertical orientation of the lander it 
resulted in dynamic controller. An analytical guidance 
law is proposed by approximating the cost function using 
the trigonometric series to address the vertical landing 
constraint (Afshari et al., 2009). The basic formulation 
has been subsequently revised to account for the crew vis- 
ibility and sensor capability (Lee, 2011). This approach 
resulted in the guidance law in which final orientation 
angle depends on the initial condition of the lander. So 


to ensure demanded orientation, a set of initial conditions 
has to be achieved, which may not be feasible since 
maneuver error results in trajectory dispersions. A sum- 
mary of all these developments around the explicit guid- 
ance law originally proposed in (DSouza, 1997) can be 
found in a good review paper (Guo et al., 2011). An aug- 
mentation is proposed in the objective function as soft 
constraint to address the requirement of lander’s terminal 
vertical orientation. Here, the aim is to minimize not only 
the energy but also the terminal error between the acceler- 
ation achieved and acceleration specified by the designer 
(Ramkiran et al., 2016). 

When guidance law derived from linearized dynamics is 
implemented for the nonlinear system, the trajectory gener- 
ated results in performance degradation. So researches are 
interested in finding the non-linear analytic guidance law 
for addressing moon landing problem. The idea is to 
approximate the position and velocity trajectory by a poly- 
nomial function of time. The coefficients of this polynomial 
are obtained using the boundary conditions. Using the 
non-linear dynamics the explicit guidance law is derived, 
which results in trajectory satisfying position, velocity, atti- 
tude boundary conditions (Sachan and Padhi, 2016). 
Though this technique enforces the boundary conditions, 
the minimum fuel consumption criterion is not addressed 
by this approach. Addressing the fuel consumption 
through the inverse polynomial guidance approach have 
been proposed that converts the landing problem into sta- 
tic optimization problem (Banerjee and Padhi, 2015). The 
coefficients of the guidance law are found by solving this 
problem using non-linear programming approach. How- 
ever, owing to the usage of simplified linearized dynamics 
and polynomial approximation of the trajectory in the pro- 
cess of mathematical derivation, the explicit solution and 
its variants, even though appear to be elegant, do not lead 
to true optimal solution in general. More importantly, path 
and hardware constraints (which are typically necessary for 
a real mission) have not been imposed and hence the solu- 
tion is not truly elegant. 

One way of including this path constraint in optimal 
control formulation is by using direct trajectory optimiza- 
tion method. Once the guidance history satisfying all the 
mission constraints are obtained in ground, it is stored 
on onboard and followed. In literature, lunar landing prob- 
lem is solved by control parameterization technique and 
time scaling transform (Zhou et al., 2010), which results 
in nonlinear guidance law addressing soft landing and ter- 
minal attitude constraint. Also, the effect of choosing the 
de-orbiting altitude on the landing site selection have been 
reported (Park and Tahk, 2011). Though the proposed 
method addressed path constraints it is not suited for solv- 
ing the multi phase lunar landing problem. The idea of 
splitting the trajectory into de-orbit, braking and vertical 
descent for addressing terminal landing accuracy and solv- 
ing this problem using Legendre pseudo-spectral technique 
(discussed in detail in Appendix A) has been handled ear- 
lier (Hawkins et al., 2006). However, mission constraints 
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resulting from landing site selection and camera look angle 
have not been addressed in that study. The first and second 
authors of this paper along with different co-worker have 
solved the powered descent phase of moon landing prob- 
lem using Legendre pseudo-spectral and Newton method 
(Mathavaraj et al., 2012). This paper compared the advan- 
tages of direct and indirect optimization methods. How- 
ever, the landing problem was solved as a ‘single phase 
problem’ considering only a limited number of mission 
constraints. 

Major technological objectives of a realistic soft landing 
in a lunar mission however are (i) to assure precise landing, 
(11) a soft touch down and (iii) selection of a favorable tar- 
get site for landing (with minimum hazard). The inclusion 
of these objectives divides the single powered descent phase 
into a multi-phase problem. For example, the ‘attitude 
hold phase’ is made necessary to ensure the appropriate 
landing site arrival and if not, course correct to reach the 
landing site suitably. This phase has necessitated that guid- 
ance command generating algorithm should also account 
for proper attitude look angle for the on-board vision cam- 
era. However, image processing of on-board camera 
requires finite processing time, which demands an 
attitude-hold phase. This in turn imposes a constraint on 
the time that must be attributed to the attitude hold phase 
so that the image processing algorithm can give a meaning- 
ful solution. All these requirements demand a multi-phase 
formulation. Inclusion of this additional constraint along 
with soft landing constraint, the powered descent phase is 
split into three phases and solved in a two-dimensional 
framework and presented in an international conference 
(Mathavaraj et al., 2016). Following similar philosophy, 
in this paper the solution is upgraded to the three- 
dimensional framework along with necessary mission con- 
straints in each phase (which include additional constraints 
to ensure zero velocity and attitude deviations in cross 
plane). Monte-Carlo analysis of optimal solution has been 
carried out to get an estimate of error ellipsoid. Robustness 
of the Legendre Pseudospectral technique is validated by 
envisioning to reach the same safe landing site for all the 
lander initial conditions taken from this error ellipsoid. It 
has been shown that the pseudo spectral method is quite 
efficient (Ross and Karpenko, 2012) and even in a multi- 
phase trajectory setting, the spectral method is found to 
be efficient as well and satisfies the terminal conditions 
accurately. 

The rest of the paper is organized as follows. In Sec- 
tion 2, a point mass based governing equations in spherical 
coordinates and a detailed formulation of converting 
multi-phase problem into non-linear programming prob- 
lem have been presented. This non-linear programming 
problem is solved for nominal condition and using this 
solution, error ellipsoid study is carried out for dispersed 
initial conditions. The guidance history in ‘braking with 
precise navigation phase’ is generated again to study the 
effect of these dispersion on fuel consumption. Simulation 


results and the mission constraints selected have been pre- 
sented in detail in Section 3. Finally, conclusions of the 
work is consolidated and presented in Section 4. In Appen- 
dix A, a generic theory of Legendre pseudo-spectral 
method is provided. 


2. Multi-phase problem formulation for lunar soft landing 


To initiate descent trajectory, first the spacecraft is to be 
de-boosted from the circular parking orbit of 100 km to 
18 km altitude (which is selected based on lunar terrain 
considerations) using a Hohmann transfer orbit. At the 
perilune height of 18 km, the lander is to enter into the 
powered descent phase. Near this phase, the propulsive 
engines are turned on and the lander velocity is decreased 
in a controlled way, to enable it to soft land on the lunar 
surface. 


2.1. Phases in moon landing 


There are three distinct phases evolved that need to be 
realized in order to land the lunar module on to the surface 
of moon. They are: (1) De-orbit maneuver phase (ii) Trans- 
fer orbit phase or coasting phase and (iii) Powered descent 
phase. The parking orbit of the lunar spacecraft is selected 
as 100 km altitude for this study. Otherwise, the method 
developed is general in nature and can be applied to any 
other orbit as well. 

Deorbit phase is to maneuver the lander from 100 km 
parking orbit to reach 18 km altitude using Hohmann 
transfer principles. The lander is designed to coast to 
18 km in transfer orbit phase, with appropriate conditions 
for the lander to approach landing site when decelerated. 
This is carefully planned such that retargeting and out- 
of-plane maneuver that costs additional fuel consumption 
is avoided. 

Powered descent phase is initiated from the terminal 
conditions of the coasting phase. This phase is further split 
into three phases to account for various mission constraints 
like camera look angle, image processing time and re- 
targeting to a safe landing site as shown in Fig. |. In this 
paper, powered descent phase is analyzed assuming that 
both deorbit maneuver phase and transfer orbit phase is 
carried out such that initial conditions for powered descent 
phase is planned as desired. 


2.2. System dynamics 


In order to describe the motion of the lander it is neces- 
sary to define a suitable coordinate system in moon cen- 
tered inertial frame XYZ as shown in Fig. 2a and 
formulate equations of motion in accordance with physical 
laws governing the system. Considering the vehicle as a 
point mass within the Moon’s gravity field, the three 
dimensional equations of motion with rotational velocity 
of moon is given by 
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Fig. 1. Powered descent phase mission scenario. 


(a) Spacecraft orbiting around the moon 
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(b) Thrust vector in body frame 


Fig. 2. Co-ordinate frames. 
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In Eq. (1), r represents the radial distance from the cen- 
ter of the moon and 0, ¢ represents the longitude and lati- 
tude traced by the lander, respectively. u,v, w represents the 


tangential velocity, across velocity and vertical velocity 
components, respectively and m represents the instanta- 
neous mass of the lander. Also pw is the moon’s gravita- 
tional constant and /,, is specific impulse of the engine 
considered for simulation. The control parameters consid- 
ered are thrust 7, thrust vector declination angle f and 
thrust vector right ascension angle « in body frame xyz is 
shown in Fig. 2b. For further details one can refer to 
Vinh et al. (1980). 

For the optimization algorithm, normalization is 
required for numerical conditioning of the variables. So, 
using basic variables like lander’s initial distance (r;) from 
the center of moon, radius of the moon (r,,) and lander’s 
initial mass (m;), normalizing variables are chosen as given 
below: 


tion. Adv. Space Res. (2017), https://doi.org/10.1016/j.asr.2017.09.016 


Please cite this article in press as: Mathavaraj, S., et al. Constrained optimal multi-phase lunar landing trajectory with minimum fuel consump- 


S. Mathavaraj et al./ Advances in Space Research xxx (2017) xxx-xxx 5 


2.3. System dynamics as algebraic constraints 
The dynamical equation constraint in each phase is for- 


mulated into non-linear programming constraints as 
explained in Appendix A 
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2.4. Cost function selection 


As discussed earlier, the spacecraft fuel consumption is 
to be minimized if it is to carry a enhanced payload mass 
for useful lunar explorations. So the cost function (J) for 
the optimal control problem is selected as minimization 
of acceleration. 


i T/m dt (3) 


2.5. Path terminal constraints 


As the trajectory is segmented into a multi-phase trajec- 
tory, there are several path constraints that are to be satis- 
fied at initial/terminal condition of the segments. Fig. | 
provides the segmentation and mission scenario of the tra- 
jectory design. 


2.5.1. Braking with rough navigation 

This phase starts from the terminal conditions of the 
coasting phase. The objective of this phase is to decrease 
the orbital velocity and to ensure terminal conditions with 
appropriate altitude and attitude for camera imaging. 
Camera imaging capability determines the terminal altitude 
of this phase. Since image resolution improves as altitude 
decreases, ground test specifies an altitude upper limit 
below which image clarity is better for identifying a safe 


landing site. Also, the lander’s terminal attitude of this 
phase is specified based on the mounting angle of the lan- 
der camera and its field of view. 


ri(tii) —Pis, Wiltz) —Wis, Ui(tis) — Wis, Vi(ti) — Vis)|” =0 
T 
ln (ty) —Fisr By (tir) - Pisy %] (tir) ~ 2s,| =v 
[1 (ty) Wis, Uy (ty) Us, VY (tip) - Y15,| : <0 (4) 


In Eq. (4) first two equation enforces the difference in 
initial (71 (t1;), Wi (tr), U1 (tii), U1 (t1)) and final 
(His Bis 2 Hsp)» states and controls to zero with respect to 
desired nominal condition (715,, Wis; Mis; Uts;s/Isp5 Bises Oise). 
The last equation enforces that the terminal velocity 
achieved should be less than the velocity sensor maximum 
capability (Wis,, Misys Uis,)+ 


2.5.2. Attitude hold 

The objective in this phase is to image the selected por- 
tion of the lunar surface, process it for the current location 
information so that in subsequent phase the guidance com- 
mand guides the lander to the targeted safe landing site. 
The ‘braking with rough navigation phase’ specifies the ini- 
tial conditions of this phase. It is to be noted that, at the 
terminal condition of ‘braking with rough navigation 
phase’, the lander is forced into a favorable orientation 
for imaging by the camera. Throughout this phase, the 
thrust firing and attitude are held constant to avoid chat- 
tering during image capturing through the following 
constraints. 


[T2(t) — Ts (tyr) Bo(t) — Bi (tip) 900) — on (n)]" = 0 (5) 


Here, the thrust vector throughout this phase 
T2(t), Bo(t), %2(¢t) remains same as the phase 1 terminal 
thrust vector T),,, B Isp? Usp SO that continuity of the control 


vector is enforced as constraints. 


2.5.3. Braking with precise navigation 

At the terminal part of the ‘attitude hold phase’, the cor- 
rect update of the position of the lander with respect to safe 
site is known. The objective of the ‘braking with precise 
navigation’ is to guide the lander from the terminal condi- 
tion of the ‘attitude hold phase’ to an altitude of 100m 
above the landing site with zero kinetic energy. Further, 
in the terminal condition, the lander’s orientation should 
be vertical, which favors the lander’s leg to touch the 
moon’s surface satisfying the constraints given below. 


[r3 (ts) —Trx 93 (t3,) — 03, 3 (ts) = os]! =0 
[ u3 (ts) —U3s U3 (ts) — 35 W3 (ts) = Way ]° =0 (6) 
EE 

[ Bs (tsp) — Bs. %3 (tar) — 035] =0 

Here, the final latitude and longitude constraint enforces 
the difference between the phase 3 final position 03(t3,), 
3(ts) and identified safe site position 03,, 3, as zero i.e. 
the lander’s terminal position is above the safe site. 
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3. Results from simulation experiments 


Extensive simulation experiments have been carried out 
for solving Multi-phase lunar landing problem explained in 
the earlier sections. The algorithm uses the suite of mathe- 
matical nonlinear programming DIDO solvers in 
MATLAB (Ross, 2007). The solver uses sequential quadra- 
tic programming which solves a sequence of subproblems 
based on the quadratic model of the original problem 
(Ross and Fahroo, 2004). The Lagrangian is constructed 
using the quadratic model with linear and non-linear con- 
straints using Karush-Kuhn-Tucker necessary condition. 
The idea behind this approach is to model at the current 
iterated solution by quadratic programming subproblem, 
then use the minimizer of this subproblem to define a next 
iterated solution. The step size at every iteration is based on 
the lagrangian merit function with a linesearch method 
(Nocedal and Wright, 2006). The difference between the 
current and previous iterated solution is checked for the 
tolerance specified. If the difference is less than the toler- 
ances, the iteration stops otherwise the search continues 
in the descent direction calculated by sequential quadratic 
programming. The covector mapping principle eliminates 
the curse of sensitivity associated with solving costates of 
the optimal control problem. The optimality of the gener- 
ated spectrally accurate solution is verified by Pontryagin’s 
Minimum Principle. 

Here, the problem is formulated as detailed in Section 2 
which consists of a set of equality and non-equality 
bounds. Normalization is performed so that the algorithm 
searches for feasible solution in the same non-dimensional 
domain as discussed above. 

For the analysis, the state and control guess history for 
each phase is simulated using the explicit analytical guid- 
ance solution derived from the linearized dynamics. Here 
the minimization of jerk is taken as cost function to obtain 
a new optimal guidance law using calculus of variation 
approach. Note the algorithm not only enforces the final 
position and velocity but also acceleration at final time. 
This gives the flexibility to ensure the desired thrust angle 
orientation at the terminal condition of each phase. In 
explicit guidance scheme, the final time is derived from 
the sixth order time-to-go polynomial which is a function 
of states (Mathavaraj and Padhi, 2017). 

However, the optimum control profile from the solver is 
used to propagate the system dynamics with suitable time 
step by Runge-Kutta fourth order method for validating 
the trajectory. In order to study the additional fuel require- 
ment due to orbit determination error as well as maneuver 
error, the error ellipse by Monte-carlo simulation is also 
carried out using dispersed initial conditions. 


3.1. Optimal control solution for nominal case 
Initial and final condition for the lander is provided in 


Table 1. Table 2 gives the parameters used for the optimal 
control problem simulation. In the formulation, these 


Table | 
Terminal Condition (T.C.) for multi-phase landing problem. 
Phase Variables Initial Final 
Vertical velocity (v1;) 0 Free 
Tangential velocity (w1;) 1.69 km/s Free 
Across velocity (w1;) 0 km/s Free 
1 Altitude (715 — rm) 18.2 km 7km 
Thrust right ascension angle («1;) Free 0 deg 
Thrust declination angle (f,) Free 140 deg 
Vertical velocity (v2;) Phase 1| T.C. Free 
Tangential velocity (u2;) Phase | T.C. Free 
Across velocity (w2;) Phase | T.C. Free 
2: Altitude (725 — rm) Phase 1| T.C. Free 
Thrust (72;) Phase 1 T.C. Vt 
Thrust right ascension angle (a2,) 0 deg Vt 
Thrust declination angle (/>,) 140 deg Vt 
Vertical velocity (v3;) Phase 2 T.C. 0 
Tangential velocity (u3;) Phase 2 T.C. 0 
Across velocity (w3;) Phase 2 T.C. 0 
3 Altitude (73. — rm) Phase 2 T.C. 100 m 
Downrange (703s) Phase 2 T.C. 10000 m 
Crossrange (1m3,) Phase 2 T.C. Om 
Thrust right ascension angle («;3;) Phase 2 T.C. 0 deg 
Thrust declination angle (/3,) Phase 2 T.C. 90 deg 
Table 2 
Mission parameters for simulation. 
Variables Values 
Gravity 1.62 m/s” 
Gravitational parameter 4902.779 x 10° km?/s? 
Moon radius 1737.1 km 
Specific impulse 320s 


conditions are absorbed as constraints in the optimal con- 
trol algorithm. As seen in Table 1, each phase has its own 
mission constraints. Phase | (‘braking with rough naviga- 
tion’) terminal condition in altitude and attitude has been 
selected based on camera sensor as explained in Sec- 
tion 2.5.1. Phase 2 (‘attitude hold’) constraint emphasizes 
on avoiding chattering during imaging. Phase 3 (‘braking 
with precise navigation’) constraints enforce soft vertical 
landing on desired site. Detailed description of these con- 
straints has been discussed in Section 2. Four main 800 N 
engine is mounted diagonally opposite manner which pro- 
vides total thrust of 3200 N. In phase 1 all four engines are 
ON to meet the thrust requirement and in phase 2,3 two 
diagonally opposite thruster are switched OFF since lan- 
der’s mass reduces due to fuel consumption in phase 1 in 
turn result in thruster requirement reduction in phase 
2,3. Table 3 provides the details of the thrust control con- 
straint in each one of the phases and should be suitably 
handled in the trajectory as a constraint. 


Table 3 
Thrust constraint for multi-phase landing problem. 
Phase Upper bound Lower bound 
1 100% of 800 x 4N 40% of 800 x 4N 
2,3 100% of 800 x 2N 40% of 800 x 2N 
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Fig. 3a shows the decreasing altitude profile of the lan- 
der satisfying all mission constraints i.e. first phase terminal 
constraint of 7 km and third phase terminal constraint of 
100 m. Fig. 3a also shows the down range traced by the 
lander and Fig. 3b shows the cross range traced by the lan- 
der reaching desired landing site at final time. Fig. 3c-e 
shows the velocity components of the lander satisfying 
the mission constraints of soft landing. It is to be noted 
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(e) Across velocity profile of the lander 


that terminal velocity attained at the ‘braking with rough 
navigation phase’ is optimal since it is not constrained in 
the formulation. To demonstrate this fact, simulation exer- 
cises have been carried out to compare fuel consumption 
with obtained optimal rough braking terminal states to 
non-optimal states. 

Case 1 and case 2 are selected such that the total 
achieved velocity is more and less respectively than optimal 
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d) Radial velocity profile of the lander 
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Fig. 3. Nominal case solution. 
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Braking with rough navigation phase terminal states. 


Case study Component Velocity Total propellant consumed 
Radial —53 m/s 
Optimal Along track 259 m/s 428.21 kg 
Across track —0.5 m/s 
Radial —15 m/s 
Case 1 Along track 300 m/s 438.56 kg 
Across track 20 m/s 
Radial —55 m/s 
Case 2 Along track 230 m/s 432.89 kg 
Across track —20 m/s 


values required to cover the desired downrange value. 
Table 4 shows the fuel savings of 10 kg if the states are free 
to take the optimum values that extremizes the perfor- 
mance index. The first simulation run depicts a case where 
the total velocity is more and the excess kinetic energy has 
to be reduced before reaching the landing site which 
demands more fuel consumption. On the other hand, the 
second simulation run depicts a case where the total veloc- 
ity is low and the lander needs extra kinetic energy for trav- 
eling (extra propellant consumption) to reach the landing 
site. It is to be noted that the multi-phase problem can be 
solved either by considering each phase separately or by 
handling them together. When solving individually, the 
process becomes very laborious and tedious especially 
when many deboost initial conditions are to be studied. 
However, in the formulation presented in this study, the 
process is carried out with all three phases together which 
makes it elegant. Fig. 3f provides the demanded thrust 
throttling percentage after satisfying the mission con- 
straints as given in Table 3. In Attitude hold phase, since 
the throttling percentage and orientation is held constant, 
the vertical velocity components increases as seen in 
Fig. 3d. In the subsequent braking with precise navigation 
phase, the throttling percentage and orientation is derived 
from the current states to achieve soft landing at landing 
site. 
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(a) Thrust Declination Angle profile ensur- 
ing vertical landing of lander 


The thrust vector declination angle profile provided in 
Fig. 4a satisfying constraints viz attaining terminal con- 
straint of 140 deg at the ‘braking with rough navigation 
phase’ and 90 deg at the ‘braking with precise navigation 
phase’ ensures the vertical landing over the desired site. 
Fig. 4b shows the thrust vector right ascension angle satis- 
fying terminal constraints. Note that in ‘attitude hold 
phase’, the throttling percentage and thrust declination as 
well as right ascension angle is maintained constant 
throughout as per the design and which is demanded by 
Eq. (5). 

In Legendre Pseudospectral method since all dynamic 
and terminal constraints are converted into either equality 
and or inequality constraints, the algorithm is capable of 
accounting additional constraints due to desired landing 
site as well. In the entire simulation, the total number of 
Legendre-Gauss-Lobatto points chosen are 70 out of which 
30 each has been allocated for ‘braking with rough naviga- 
tion’ and ‘braking with precise navigation’ phases and the 
remaining 10 is used for ‘attitude hold phase’. Since the 
acceleration is held constant in ‘attitude hold phase’ the 
number of grid points is chosen less when compared to 
other phases (Ross and Fahroo, 2004). The numbers of 
grid points are sufficient since the control profile is used 
to propagate the system dynamics and the trajectory has 
been validated. Hence it is shown that using pseudo- 
spectral knots, multiphase lunar landing problem is solved 
efficiently satisfying all mission constraints with a single 
optimal control framework. Table 5 signifies the reduction 
of fuel consumption due to multiphase formulation when 
compared to explicit guidance (Mathavaraj and Padhi, 
2017) for soft landing at desired site. 


3.2. Study on effect of error dispersions 

In this design, from ‘deboost phase’ to the end of ‘rough 
braking navigation phase’, the navigation process is 
planned to rely on inertial navigation system totally. The 
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Fig. 4. Nominal case solution. 
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Table 5 
Comparison with explicit guidance. 


Case study Initial position Values Propellant consumed 
Altitude 18 km 
Legendre Downrange —389 km 430.91 kg 
Pseudospectral Cross range —600 m 
Altitude 18 km 
Explicit Downrange —389 km 454.66 kg 
Cross range —600 m 


from laser altimeter, position reference from hazard detec- 
tion camera and velocity components from velocity meter 
are obtained from sensors that are switched on only at 
the start of the ‘attitude hold phase’. The sensors are so 
designed that they can perform efficiently only in this alti- 
tude range. At the end of 35s period (or start of ‘precision 
braking navigation phase’), it is expected that the absolute 
value of position, altitude and velocity can be estimated by 
the onboard process. The error estimates and the capability 
of the algorithm to counter the error for precision landing 
are explained in the following. When deboost of the lander- 
craft from the 100 km circular parking orbit is initiated, it 
is effected based on the Orbit Determination (OD) results 
for that parking orbit. The errors of OD estimation along 
with other factors such as maneuver performance and 
misalignment of thrusters have a considerable influence 
on the precision of descent trajectory of the lander and 
therefore, a simulation study is conducted using Monte- 
Carlo procedure. Based on Chandrayaan-1 on-orbit experi- 
ence, the OD estimation statistics are determined to be 
1 km (3¢) in position and 20 cm/s (3a) in velocity and 
for a 100 km circular parking orbit. The corresponding 
36 errors in radial, along and across track directions is pro- 
vided in Table 6. Using the statistical values of Table 6, one 
hundred initial conditions of OD for 100 km parking orbit 
are randomly generated for deboost conditions and the de- 
boost is carried out. 

The de-boost orbit is selected as 100 x 18 km elliptical 
transfer orbit based on the mean Moon model and the 
expected landing site undulations. The deboost dV maneu- 
ver error is assumed to be 1% of total thrust magnitude 


4x 800N and about 0.1 deg in thrust direction. A one 
thousand different sets of state vectors accounting for 
thrust and alignment errors are generated from the one 
hundred initial conditions generated earlier. The trajectory 
is propagated for the nominal coasting duration in transfer 
orbit phase to reach 18 km perilune position and the result- 
ing deviations in altitude, down and cross ranges from such 
initial state vectors are generated. When these initial condi- 
tions are used in open loop simulation during ‘rough brak- 
ing navigation phase’, the terminal errors at the end of 
‘attitude hold phase’ are accumulated to be 2 km in alti- 
tude, 4km in downrange and the 0.1 km in cross range. 
The statistical details are provided in Table 6. 


3.3. Regeneration of braking with precise navigation phase 


At the end of the ‘attitude hold phase’, due to the esti- 
mation process onboard for altitude, position and velocity, 
there is a fairly good estimate of the errors in the initial 
conditions available for correction. Using these initial con- 
ditions and knowing the final landing site details, the opti- 
mal trajectory are recomputed using Legendre pseudo 
spectral method for all one thousand cases and all states 
of the simulation are plotted in Fig. 5. Fig. 5a provides 
the dispersions of altitude, downrange and crossrange aris- 
ing due to various error sources considered above. It is to 
be noted that the error sources considered have a signifi- 
cant effect on the final absolute velocity of the lander. 
Fig. 5b portrays that in spite of all the velocity error noted 
in such cases, the re-computation of the trajectory is 
enacted accurately such that the final velocity is almost 
zero for all cases satisfying soft landing mission constraint 
but at the cost of some extra fuel spending. However, based 
on the initial position and velocity with respect to the land- 
ing site, the acceleration demand on the lander to reach the 
landing site is found to vary. Table 7 provides the addi- 
tional fuel requirement due to dispersion in the states but 
managed the fuel expenditure within the fuel margin. 

Fig. 5c and d gives the percentage of the thrust 
demanded and fuel mass consumption for all the cases con- 
sidered. Note that the demanded thrust requirement is 
within the main engine capability of the lander of 


Table 6 
Orbit determination error. 
Phase (Altitude) Component Position Velocity 
De-orbit (100 km) Radial 400 m 6 cm/s 
Along track 900 m 18 cm/s 
Across track 10m 2 cm/s 
Braking with rough navigation (18 km) Radial 0.5 km 3 m/s 
Along track 3km 0.3 m/s 
Across track 0.05 km 0.09 m/s 
Braking with precise navigation (5.5 km) Radial 2km 4 m/s 
Along track 4km 1 m/s 
Across track 0.1 km 0.2 m/s 
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Fig. 5. Regeneration of braking with precise navigation phase. 


2 x 800 N. Fig. Se shows the demanded thrust vector decli- 
nation angle profile whose initial condition is same as ‘atti- 
tude hold phase’ angle 140 deg satisfying the continuity at 
the transition from ‘attitude hold phase’ to ‘braking with 
precise navigation phase’. The terminal condition of the 
thrust declination angle is 90 deg for all error cases ensur- 


ing vertical orientation at the final time. Fig. 5f shows the 
thrust vector right ascension angle profile demanding the 
corrective out of plane maneuvers wherever the crossrange 
error has to be nullified. It is to be noted that at the final 
time, the thrust vector right ascension angle is 0 deg which 
favors the lander’s leg orientation for vertical touchdown. 
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Table 7 
Fuel consumption in braking with precise navigation phase. 


Case Error component Position Fuel consumption 
Radial 400 m 
Nominal Along track 900 m 60.47 kg 
Across track 10m 
Radial 2km 
Case 1 Along track 4km 65.35 kg 
Across track 0.1 km 
Radial —2km 
Case 2 Along track —4km 68.58 kg 
Across track —0.1 km 


4. Conclusion 


A unified multi-phase constrained optimal trajectory 
problem is posed and solved to optimally guide a lander 
from the perilune of the transfer orbit to a desired altitude 
for smooth vertical descent in order to achieve soft-landing 
on the surface of the moon. This design accounts for vari- 
ous operational constraints and suitably segments the tra- 
jectory into three phases namely, the ‘braking with rough 
navigation phase’, ‘attitude hold phase’ and ‘braking with 
precise navigation phase’ of the landing sequence. All rele- 
vant mission constraints such as upper and lower bounds 
of the thrusters, attitude hold constraint to cater for the 
look angle constraint of the downward looking vision cam- 
era and time required for image processing have been con- 
sidered in the optimal control problem formulation. 
During the cases of high terminal error due to dispersions, 
the ‘precision navigation terminal phase’ guidance history 
is regenerated and the additional fuel requirement is 
revealed. The proposed approach also leads to a minimum 
fuel consumption solution, which is an important objective 
of this study. Extensive simulation studies show that fol- 
lowing the strategy proposed in this paper the desired ter- 
minal position, velocity and attitude are achieved. In 
addition, proper orientation and safe soft landing is 
ensured over the desired landing site. It is also shown that 
the dispersions due to orbit determination and maneuver 
are absorbed and a precise landing can be achieved. 
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Appendix A. Legendre pseudospectral method: a brief 
summary 


Even though the Legendre pseudo-spectral method has 
evolved as a popular technique and is used by many 


researchers across the globe, the theoretical details of the 
method is not found in most of the literature to the best 
of the knowledge of the authors. Hence an attempt has 
been made here to explain the technique in a lucid manner 
so as to make the technique accessible and understandable 
to large practicing engineering researchers and to inform 
what precautions one must adhere to operating successfully 
while solving practical problems. 


A.1. Problem objective 


A generalized optimal control setup is with the objective 
to minimize the cost function J in Eq. (A.1) 


J=o(x(t).t) + [LW U0,dd 


i) 


(A.1) 


subjected to the following dynamical constraints in Eq. 
(A.2) and trajectory constraints Eqs. (A.3), (A.4) 


X(1) = f(X(), U(9) (A.2) 
e(X (to), X(tp)) = 0 (A.3) 
h(X(t), U(t)) <0 (A.4) 


where f,¢; represents the initial and final time, X repre- 
sents the state vector, U represents the control vector, 
o(X(ty), ty) represents the terminal cost term, L(X, U, t) 
represents the integral cost term, A(X, U), e(X(t)X(t,)) 
represent the path and terminal constraints, respectively. 


A.2. Discretization of the system dynamics 


Usually collocation methods use piecewise-continuous 
functions to approximate the state and control variables 
at the arbitrary subintervals. However, the Legendre 
pseudo-spectral method uses globally interpolating 
Lagrange polynomials as trial functions to approximate 
the state and control variables at the nodes, which are 
called Legendre-Gauss-Lobatto (LGL) points. 


A.2.1. Selection of grid points 
The time interval (ty—¢;) can be transformed to t 
domain by change of variables 


ti)/ (ty — ti) 


This transformation from the real time domain to the t 
domain is required because the LGL points are defined 
between [—11]. So the N+1 points are t;/=0,...,N 
which specifies t) = —1,ty = 1 and other LGL points are 
the zeros of the derivative of the Legendre polynomial of 
degree N, (Ly). 


tO 14 2(¢ (A.5) 


A.2.2. State dynamics approximation 

As stated earlier, approximating the state X(t) and con- 
trol U(t) trajectory by N” degree Legendre polynomial 
leads to 
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X(t) A") = SOX (Dl) 
U(a) © UN) = DUC) i) (A.6) 


where, Lagrange polynomials p,(t) in t domain is given by 
g(t) 

t) => > Fr A.7 
pil ) (t _ 11) 2(t1) ( ) 
where g(t) define the LGL grid points and it’s correspond- 
ing derivative is g(t) 

g(t) = (1— 1) Ly(t) 


g(t) = —N(N + 1)Ly(t) (A.8) 


Eq. (A.6) represents polynomials of N” order, which is 
used to interpolate the functions at the LGL points. Substi- 
tuting Eq. (A.8) in Eq. (A.7) and simplifying provides 


1 (1—7?)Ly(t) 
pit) = N(N + 1)Ly(t) * TOT a) 
1 if l=k 
nile) = {4 if 1#k 


Differentiating Eq. (A.6) once and substituting in Eq. 
(A.2) provides 


N 
X(t) = SOX (41) pi(te) (A.10) 
0 
aa. a 1 (1 — ?)Ly(z) 
a ay dt | N(N + 1)Ly(t1) * TT elt) 
For t = %;k 4 1 Eq. (A.11) is rewritten as 
1G (A.12) 


~ Lw(ti)(te — 7) 


For t = %;k = / results in an indeterminate form. Using 
L’Hopital’s rule and carrying out necessary substitutions 
result in 


. Ly'(t) 
~ A.13 
pi(t) 2Ly(t7) ( ) 
Caset=t; k=1#0,N;  pi(te) =0 
Caset=%; k=1=0; pi(t) = a 
Caset=%; k=I=N; pi(te) = NOCD 


The differentiation matrix Dy = ;(t,) in Eq. (A.14) is 
an (V+1)x(N+1) matrix after the discretization of 
the right-hand side of the state equation. 


Ly (te) 


aa ot 
N(N+41) a oe 
D= [Dal = 4 cs (A.14) 
ae == h 
0 otherwise 


A.3. Discretization of the cost function 


The idea is to convert the integral term in the cost func- 
tion into an algebraic constraints using Gauss-Lobatto 
integration rule. Using Eq. (A.6) in Eq. (A.1) and after sim- 
plifying one obtains 


JN = p(X", tw) 


1 


4 Gedy cn), U(t:), Tt) / [pi(t)|dt 


1 


(A.15) 
Let the weight factor be 
1 
m= f p,(t)dt 
=f 
Lf (ga) ) 
= dt A.16 
eo aie 


In order to evaluate w;, consider Christoffel Darboux 
Identity as given below. 


Bn(t,Z) = y a 


k=0 VK 


_ [Senet a ee) 
(Amt? n(T ~~ z)/Am) 


Evaluate Eq. (A.17) for N,N —1 degree polynomials 
and using Eq. (A.7) in Eq. (A.17) results in 


N+1 


(A.17) 


By(t,z) + By-1(t,Z) 
_ QN +1) _ Zy(z)g(t) — Lw(t)g(z)) 
iam, ace <= (A.18) 


Substituting z = t ie. zeros of g(t) in Eq. (A.18) and 
integrating the equation on both sides once results 


1 
N+1 
/ (Bu(c.2) + V By-a(tst))d 
-1 


_QN+) ' g(t) 
ON ints) | ESS 


(A.19) 


Evaluating right-hand side of Eq. (A.19) provides 
2N +1 : 
= ove) x Ly (t.)g! (th) We 
(2N +1)(N +1) y 
2 


Similarly, evaluating left-hand side of Eq. (A.19) results in 


: N+1 
/ (But Th.) + ai 
=f 


NEL 
+O ff Bualendde 


[Lv (te) ]’we (A.20) 


1 
By-a(svn) a= ff By (t, %%)at 


1 


(A.21) 


Using the identities 


Lesa(t) Lea] =0 


k>0 
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[Qn + DEn(@) = tyes) = bya) 


evaluating first part of Eq. (A.21) with Eqs. (A.7), (A.17) 
results in 


/ By(t, |) dt = ‘oC [—Li(t) + L1(1)]|, 


(A.22) 
Next, evaluating Eq. (A.21) using Eq. (A.22) provides 


1 1 
/ By(t,z)dt + (=) i By-1(t,z)dt 
= N “| 


_ (2N +1) 
Sr a (A.23) 
Using Eqs. (A.20), (A.23) in Eq. (A.19) results 
2N+1)(N+1 2N +1 
ON ND x thy (ea)PPme = “AY 
from which w,; is evaluated as 


N(N + D[Ew(t)) 


So the optimal control problem represented by Eqs. (A.1)— 
(A.4) is approximated by the following nonlinear program- 
ming problem. The performance index is reformulated as 


ty — to) & 
JY = p(X" ty) + We) 14x(a), U (tk), te) We 
k=0 
2 1 
N(N +1) [Ly (te)? 
the state equation as 


(ty — to) 
2 


and the constraint equations as 
e(X(t),X(tw)) = 0 
h(X(t,), U(te)) < 0 


Wk 


I (X(t), U(t), te) — Ck = 0 (A.25) 


Legendre Pseudospectral method is mathematically 
complex and the number of nodes used for simulation 
determines the dimensional complexity of the problem. In 
general, Legendre Pseudospectral method needs more 
nodes for accurate solution for a multi-phase problem. 
However, Legendre Pseudospectral knotting allows calcu- 
lation of solutions efficiently even with fewer nodes, which 
reduces the computational time. Here, the entire time 
domain is divided into smaller subdomains where Legendre 
Pseudospectral knots is enforcing the continuity of the 
state and control variables in the transition regions. For 
more details refer to Ross and Fahroo (2004). 
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